UHIVERSITY OF HUBOIS 
DIGITAL COMPUTER 

LIBRARY ROUTIH E M19 ^ 232 
By Gene H. Golub 

TITLE Solution of the Matrix Equation Ax * X Bx Where A and B are 

Symmetric and B Is Positive Definite. 
TYPE Entire Program. 

ACCURACY Depends on condition of B. 

DURATION (a) 33 seconds to input routine 

■z o 

(b) (101 + 2.1) n + 5.2 n milliseconds where n « order of 
matrix and I = number of iterations. I varies from k 
for n = 3 to 7 for n = 19. 
CAPACITY Matrices of order 19 . 

METHOD OF USE (l) Read program tape into memory in usual way. If no reading 

i error has been committed machine will stop on 2k (kkf)^^ 

(2) Place data tape in reader and restart machine with black 
'- switch. 

(3) After computation has been completed, results will be 
punched out for printing. 

(k) The machine stops on 2k (0^7),/- . A new problem can be 
begun by repeating step 2. 
FUNCHtJIG OF THE DATA To compute the eigenvalues or eigenvectors of the matrix equation 

Ax. = X.Bx. (i = l,2,...,n) where A and B are symmetric and B 
is positive definite proceed as follows: 

(a) Scale the matrices so that 

S 2 n ? 

Z a < 1/2 and Z b, . < l/2 

i,J-a 1J i^J-d. 1J 

(b) Punch the lower half of the matrix A, 

a n 

^1 a 22 

a 31 a 52 a 33 

a nl a n2 a nn 
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rov by row, as a sign followed by up to twelve decimal 

digits. The last element a is followed by N. 

nn " 

(c) Punch the matrix B in a similar fashion but follow the 

last element b by an N.J.F. or L. 
nn 

Let By i = 8^ i = 1,2,..., n. 

Then the results of punching these characters is as follows: 



N 
J 
F 
L 



\ . , (i = 1, ...,n) are punched 
A.., x (i = l,...,n) are punched 
5., \. (i = l,o.o,n) are punched 



5., \ , x (i=l,o..,n) are punched. 

(d) The last character is followed by a sexadecimal character 
p which determines the number of decimal digits to be 
printed. The character p can assume the values 1,2,3, ».. 
9,K,S,N where K = 10, S = 11, and N = 12. 
THE PRINTED RESULTS First, the eigenvalues of B are punched if desired. Next if 

only the eigenvalues of the system Ax = Xx. are computed then 
they are punched in a single column. If the eigenvalues and 
eigenvectors are computed, then following each X. will be the 
corresponding vector x. . 

Throughout this program, there are a number of checks . They 
are as follows: 

2 
FF 02S Matrix A fails to have £-±2: elements 

2 
FF 02N Matrix B fails to have 2_+£ elements 

FF 02J A and B fail to have same number of elements 

FF 02F At least one eigenvalue of B < 5xl0" 10 

DISCUSSION OF METHOD USED 

To solve the equation |A - XB | = we first solve |B - 8l| = 0. 

Since B is symmetric we know that an orthogonal matrix Y and a 

T 
diagonal matrix A exists such that B = YAY ; Y is the matrix 

of eigenvectors and A the matrix of eigenvalues for B. The 
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elements of A are non -negative since B is positive definite, 
so we may write 

B = YA 1 ' 2 A 1 / 2 Y 1 " = UU 31 

1/2 
where U = YA . Hence we have 

|a - xb| = |a - \mj r | =0 

and 

|lf 1 AU~ T - XI | = 0. 

-1 -T 
The matrix U AU is symmetric so that there is an orthogonal 

matrix Z and a diagonal matrix A such that 

U AU = Z/\Z . 

The eigenvalues \. are, then, the elements of A . By definition 
of an eigenvector, we have 

(U" 1 AU _T ) z, i =^ j[ z i 

T -T 
= \. IT U a. 
11 

-T T -T -T 

and AU ± z, ± = L UU U ^ = X. BU z ± 

-T 
Then, if we let x. = U fc,. we have Ax. = X.Sx. 
l i l i i 

so that \. and x. are the eigenvalues and eigenvectors of the 
system |A - XB| = 0. 

-T 
The matrix multiplication U Z is performed by replacing the 

-T -1 -T 

identity matrix by U when finding the eigenvalues of U AU • 

T 
Therefore, £. Bx. = 1 for i = l,2,.o.,n. 

NOTE Frequently, it is of interest to find the eigenvalues and 

eigenvectors of Hermitian matrices. The present routine can 
handle these problems, although in a somewhat wasetful manner. 

Suppose A = G + i H 
and B = M + i N 

Then if B is positive definite, 



- k - 

(G + i H)(x + iy) * X(M + iN)(x+iy) 
Gx .- Hy = \(Mx - Ny) 
Hx + Gy = \(Kx + tty) 



That is 



g, -ir 

B, G_ 




y 


« X. 


M, -H 

JL 21 


■ 


"x 

y 



But -H = H and 
Therefore 



-N = if 



T 
G, H 




x 


= X 


T 
M, N 




X 


2L? G_ 




_y_ 




_N M_ 


. 


UJ 



Solving this equation will yield 2n roots. Each root of the 
original system will be duplicated. 



DAT E February 13* 1957 

CODED B Y /7 j^^L — A 

APPROVED B Y?) f , pfjJiJL^ 



LOCATION 



ORDER 

00 3K 
00 F 
00 473F 
00 F 
00 F 
00 F 
00 20F 
00 F 
00 190F 
00 F 
00 253F 
00 F 
00 30F 
00 F 

00 3560F 

00 F 
00 302F 
00 F 

00 2560F 

00 F 
00 110F 

00 20K 
00 F 
00 F 
00 (p)F 
00 F 
80 F 
00 (11 )F 
00 F 
00 F 
00 F 
00 F 
80 F 
00 F 
00 (n)F 
00 (n)F 



NOTES 
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Location of parameters 

Location of P-l6 

Location of R-l 

Location of Y-l 

Location of matrix A on drum 



Location of subroutines on drum 



Location of lA-k 



LOCATION | ORDER 

00 F 
00 500J 
00 IF 
00 IF 
80 F 
00 (n 2 )F 

00 50K 

Library Routine Y-l 
(modified) 

00 71K 
50 110F 
50 L 
26 S8 
00 SS 
00 73F 
50 2L 
26 110F 
L5 S5 
80 IF 
IA 36l 
kO 27L 
L5 6S5 
k6 10L 
k6 2TL 
50 SN 
50 7L 
8 26 S8 

00 100SS 

00 165F 

00 IF 

10 50 ( )F 
50 10L 

11 26 SN 
00 IF 

12 50 110F 
50 12L 



NOTES 
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from 55 



Read input routine off drum 



Enter input routine 



Read M-4 off drum 



by 6 



Enter M-k 



Read in routine II off drum 



LOCATION 



ORDER 



NOTES 



PAGE 3 



M 19 



13 
14 

15 
16 

17 
18 

19 
20 
21 
22 

23 

2k 

25 
26 

27 

28 

29 
30 



26 S8 
00 300SS 
00 45F 
00 IF 
50 S6 
50 15L 
26 S8 
00 350SS 
oo 56F 
50 17L 
26 110F 
00 IF 
50 110F 
50 19L 
26 S8 
00 425SS 
00 70F 
50 21L 
26 110F 
00 IF 
50 SN 
50 23L 
26 S8 
00 100SS 
00 165F 
L5 37L 
kQ 8SN 
00 IF 
50 ( )F 
50 27L 
26 SN 
00 IF 
50 110F 
50 29L 
26 S8 
00 500SS 



Read in P-l6 off drum 



Enter routine II 



Read routine III from drum 



Enter routine III 



by 5,6 



Read M-k from drum 



Modify M-4 



Enter M-k 



Read routine IV from drum 



LOCATION 



31 



32 



33 



3* 



35 



36 



37 



ORDER 

00 43F 
00 IF 
50 S6 
50 32L 
26 S8 
00 350SS 
00 56F 
50 3^L 
26 110F 
24 L 
JO F 
50 27L 
L5 F 
32 130F 
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Input Re 


Dutine 




00 110K 





K5 F 




42 38L 


1 


50 S3 




50 1L 


2 


26 4lL 




L3 4f 


3 


36 4l 




FF 043F 


4 


L5 3F 




4o 8f 


5 


LO 8S5 




42 2S5 


6 


42 6S5 




00 20F 


7 


46 6S5 




50 6S5 


8 


75 6S5 




10 IF 
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Read P-l6 from drum 



Enter routine IV 



Plant link 



Input matrix A 



Compute parameters 



LOCATION 


ORDER 


NOTES PAGE 5 


9 


k2 9S5 
00 20F 






10 


k6 30L 
Lk 6S5 






11 


10 21F 








Fk 30L 






, 12 


42 4S5 
L5 30L 






13 


42 15L 
4l OF 






Ik 


F5 15L 
k2 15L 


from 19 




15 


00 IF 






16 


L5 F 
32 l6L 


by 13,1^,2: 
from 25 


'; 




kO 190S3 


by 17 




17 


F5 16L 








ko 16L 




square matrix A 


18 


F5 OF 








40 OF 






19 


LO 40L 
36 l4L 






20 


kO IF 
L7 IF 






21 


32 22L 
L5 15L 






22 


k2 30L 








L5 15L 


from 21 




23 


L4 OF 
42 15L 






2k 


L5 OF 
LO 2S5 






25 


32 15L 
F5 40L 






26 


kO 40L 
FO 2S5 
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LOCATION 
27 

28 

29 
30 
31 
32 
33 
3^ 
35 
36 
37 
38 

39 
IfO 
hi 
k2 
^3 



ORDER 

36 28L 
22 12L 
JO 190S3 
50 28L 
26 S8 
00 S9 
00 ( )F 
50 1023S3 
50 S3 
50 31L 
26 4lL 
00 38F 
40 S5 
81 4F 
00 20F 
k6 ~1S5 
L3 ^F 
32 36L 
FF O^F 
L5 8F 
LO 3F 
kO OF 
L3 F 
32 ( )F 
FF 045F 
00 F 
80 F 
00 IF 
50 42L 
26 999F 
00 F 
00 4lL 
00 F 
26 !OL 

26 IN 



NOTES 
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by 22 



Record A on drum 



Input matrix B 



by 



by 26 



LOCATION 


ORDER 




NOTES PAGE 


7 


Modified N 


-2 

00 800K 











JO 110F 
50 L 








1 


26 S8 

00 SS 




Place input routine on drum 




2 


00 73F 
L5 ^OL 








3 


kO 500F 
26 999F 
26 800N 








Routine . 


EI 
00 110K 











K5 F 
k2 kkL 




Plant link 




1 


L5 S5 
36 4L 




Print eigenvalues of B ? 




2 


92 135F 
92 515F 








3 


L5 1S5 
1+6 13L 




Plant print parameter 




k 


*U 3F 
L5 5S5 


from 1 






5 


kO kF 
L5 S3 


by 17; from 19 






6 


^0 5F 




Store eigenvalues 






40 S3 


by 16 


of B consecutively 




7 


LO 7S5 
36 9L 




Are eigenvalues of B > 5 x 10" ? 




8 


FF 046F 
2*1- OOOF 








9 


IA 7S5 
LO J4-F 


from 7 


- 




10 


32 11L 




Test for smallest eigenvalue 






L5 5F 




■ 
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LOCATION 



11 
12 

13 
Ik 

15 
16 

17 
18 

19 
20 
21 
22 

23 
24 

25 
26 
27 
28 



ORDER 



NOTES 
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40 4F 


L5 S5 


32 


15L 


L5 5F 


50 


( )F 


50 13L 


26 s6 


92 


131F 


92 


515F 


F5 


6l 


40 6l 


F5 


5L 


F4 


3F 


42 


5L 


F5 


3F 


ko 


3F 


LO 2S5 


32 


5L 


L5 


kF 


50 20L 


26 


S7 


49 3F 


19 38F 


ko 


3S5 


F5 2F 


LO 


3F 


36 27L 


L5 


3F. 


10 


IF 


ko 


3F 


F5 


3S5 


22 


22L 


41 6f 


41 


OF 


L5 


3S5 


k2 jkL 



from 10 



by 3 



from 12 



from 26 



from 24 



LOCATION 



29 
30 
31 
32 
33 
34 
35 
36 

37 
38 
39 
40 
41 
42 

43 
44 



ORDER 



41 8F 
L5 4S5 
L4 6f 

42 33L 
42 36L 
L5 S3 
32 32L 
50 32L 
26 S7 
L5 ( )F 
^0 OF 
10 ( )F 
66 2F 
S5 F 

32 36L 
40 ( )F 
L5 33L' 
L4 6S5 
42 33L 
42 36L 
F5 8F 
4o 8f 
LO 2S5 
32 33L 
F5 31L 
40 31L 
F5 6f 
4o 6f 

LO 2S5 
36 29L 
00 IF 
22 ( )F 

00 800K 
JO 110F 
50 L 
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from 43 



by 41 



by 30,38; from 40 
by 28 



by 31,38 



by 



Record routine II on drum. 



LOCATION 



ORDER 



26 S8 
00 300SS 
00 4^F 
L5 40S8 
40 501F 
26 999F 

26 800N 

Routine III 

00 110K 



NOTES 
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L5 3L 


9 


46 13L 




L5 4S5 


• 10 


L4 1L 




42 13L 


11 


41 2L 




41 ( )F 


12 


2L 13L 




S5 F 



K5 F 

42 69L 

L5 9S5 
00 20F 
46 5L 
41 OL 
50 SK 
50 3L 
26 S8 
00 S9 
00 F 
L5 13L 

46 3L 

43 14L 

47 l4L 
43 11L 
41 1L 



by 2,6 



by 2 
from 31 



from 23 



by 7,21 
from 19 



LOCATION 


ORDER 


NOTES 
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13 


50 SK 


by9 ' 17 W . T 

by 10,16 AU 






lk (u.^F 




14 


Lk ( )F 
ko ( )F 


by 7,20 
by 6,20 






15 


L5 13L 
Lk 2S5 








16 


k2 13L 
L4 8S5 








17 


46 13L 
F5 2L 








18 


40 2L 
LO 2S5 








19 


32 12L 
L5 i4L 








20 


Lk 8S5 
40 l4L 








21 


42 11L 
F5 IL 








22 


40 IL 
LO 2S5 








23 


32 8L 
4l 2L 








2k 


47 25L 
00 IF 


Waste 






25 


L5 ( )F 


by 24,27 ; from 29 






40 SK 


by 27 






26 


L5 25L 
L4 8S5 








27 


40 25L 
F5 2L 








28 


40 2L 
LO 2S5 








29 


36 25L 

F5 OL 








30 


40 OL 
LO 2S5 
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LOCATION 


ORDER 
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31 


32 5L 
F5 5S5 






32 


40 3L 

41 0L< 






33 


4l 1L 
43 39L 


from 55 




34 


47 42L 
43 42L 






35 


L5 OL 
L4 4S5 


from 50 




36 


00 20F 
46 4lL 






37 


L5 1L 
L4 40L 






38 


42 4lL 
4l 2L 






39 


00 IF 








41 (0)F 


by 33,48 




40 


2L 4lL 








S5 SK 


from 46 




4l 


50 (u ij )F 
74 (a ij )F 


by 36,44 ; 
by 38,44 


from 40 

If 1 (Atf T ) 


42 


L4 OF 


by 34,47 






40 OF 


by 34,47 




43 


L5 4lL 
L4 6S5 






44 


42 4lL 
46 4lL 






45 


F5 2L 
40 2L 






46 


LO 2S5 
32 40L 






47 


L5 42L 
L4 8S5 
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LOCATION 


ORDER 
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48 


40 42L 
42 39L 






4 9 


F5 1L 
40 1L 






50 


LO 3L 
36 35L 






51 


JO OF 
50 51L 






52 


26 s8 








oo S9 


by 54 




53 


00 IF 
F5 OL 


by 56 




54 


F4 52L 
4o 52L 






55 


L5 53L 
L4 8S5 






56 


46 53L 
F5 3L 




■ - 


57 


40 3L 
F5 OL 






58 


40 OL 
LO 2S5 






59 


36 33L 
41 L 






60 


50 S3 


by 65 






50 6ol 


from 69 




6l 


26 s8 








00 S9 


by 63 




62 


00 IF 
F5 6lL 


by 67 




63 


f4 l 
4o 6il 






64 


L5 6ol 
L4 62L 
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LOCATION 



65 



66 



67 



68 



69 



ORDER 

LO 8S5 
k6 60L 
L5 62L 
IA 8S5 
k6 62L 
F5 L 
^0 L 
LO 2S5 
32 60L 
22 ( )F 

00 800K 
JO 110F 
50 L 
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26 S8 



00 425SS 
00 70F 
L5 *K)S8 
kO 502F 
26 999F 

2.6 800N 



Routine IV 
I 
Print Routine 



00 110K 
K5 100F 
k2 38L 

L5 3S5 
\2 2L 
19 38F 
00 ( )F 
kO OF 

L5 1S5 
k6 29L 
92 131F 



by 



Record routine III on drum 



by 1,1^ 



by 15 



LOCATION 



10 
11 
12 
13 

15 
16 

17 
18 

19 
20 
21 
22 



ORDER 



50 39L 
L5 OL 
iA 29L 
kG 29L 
75 ^1L 
L5 OF 
SO F 
32 5L 
S5 F 
40 IF 
L5 OF 
50 ^OL 
66 IF 
S5 F 
32 12L 
kO kF 
L5 2L 
L4 3S5 
k2 2L 
L5 ^2L 
kO 6l 
46 4L 
F5 12L 
k2 12L 
LI 1L 
kO 1L 
36 2L 
41 7F 
92 131F 
92 515F 
50 5F 

7J S3 
54 f 

50 21L 
26 S6 
L5 S5 



NOTES 
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by 15 



00 F 00 IF 



00 F 00 10F 



00 F 00 F 



h: eigenvalue factor 
5: eigenvector factor 



from 27 



by 35 



LOCATION 



23 
24 

25 

26 

27 
28 

29 
50 
31 
32 
33 

34 

35 

36 

37 

38 

39 

4o 



ORDER 



00 IF 
32 34L 
92 151F 
L5 4S5 
L4 7F 
42 28L 
4l 8F 
00 IF 
92 131F 
92 515F 
50 4f 
7J ( )F 
54 F 
50 29L 
26 S6 
L5 28L 
L4 2S5 
42 28L 
F5 8f 
4o 8f 
LO 2S5 
36 27L 
92 135F 
F5 20L 
F4 7F 
40 20L 
F5 7F 
40 7F 
LO 2S5 
36 19L 
92 131F 
22 ( )F 
00 F 
00 IF 
00 F 
00 OF 



NOTES 
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from 33 



by 31 
by 4,6 



LOCATION 


ORDER 


NOTES PAGE 17 


kl 


00 F 
00 10F 






k2 


L4 21L 
46 21L 

00 800K 









JO 110F 
50 L 






1 


26 S8 
00 500SS 






2 


00 43F 
L5 40S8 






3 


40 503F 
26 999F 
26 800N 

00 110K 






Library Routine M-4 






Change Word 45 








26 163L 








50 6f 






Add to End 


of Code 
L3 6F 
32 110L 
L5 15F 
22 45L 

00 800K 









JO 110F 
50 L 






1 


26 S8 
00 100SS 






2 


00 165F 

L5 40S8 






3 


kO 504F 
26 999F 
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LOCATION 



ORDER 



NOTES 
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26 800N 

00 190K 

Code P-l6 

00 800K 
JO B6 
50 L 
26 S8 
00 350SS 

00 56F 

L5 ^0S8 

ko 505F 

26 999F 
26 800N 



Code 2 



00 800K 

- 7 
2k 71N 



